{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "import numpy.random as rnd\n",
    "import tsplib95\n",
    "import tsplib95.distances as distances\n",
    "\n",
    "from alns import ALNS\n",
    "from alns.accept import HillClimbing\n",
    "from alns.select import RouletteWheel\n",
    "from alns.stop import MaxRuntime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SEED = 7654"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The travelling salesman problem\n",
    "\n",
    "The travelling salesman problem (TSP) is a classic problem in operations research. It asks how to construct the minimum distance tour between a number of nodes, such that each node is visited once and the tour concludes at the starting city (that is, it forms a cycle). It is perhaps the best-known problem in the class of [NP-hard](https://en.wikipedia.org/wiki/NP-hardness) problems.\n",
    "\n",
    "## Data\n",
    "There are a considerable number of test data sets available for the TSP, varying in size from a hundred or so locations to many hundreds of thousands. For the sake of exposition, we shall use one of the smaller data sets: the data from the XQF131 VLSI instance, made available [here](http://www.math.uwaterloo.ca/tsp/vlsi/index.html#XQF131). It consists of 'only' 131 nodes, with an optimal tour length of 564."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA = tsplib95.load(\"data/xqf131.tsp\")\n",
    "CITIES = list(DATA.node_coords.keys())\n",
    "\n",
    "# Precompute the distance matrix - this saves a bunch of time evaluating moves.\n",
    "# + 1 since the cities start from one (not zero).\n",
    "COORDS = DATA.node_coords.values()\n",
    "DIST = np.empty((len(COORDS) + 1, len(COORDS) + 1))\n",
    "\n",
    "for row, coord1 in enumerate(COORDS, 1):\n",
    "    for col, coord2 in enumerate(COORDS, 1):\n",
    "        DIST[row, col] = distances.euclidean(coord1, coord2)\n",
    "\n",
    "SOLUTION = tsplib95.load(\"data/xqf131.opt.tour\")\n",
    "OPTIMAL = DATA.trace_tours(SOLUTION.tours)[0]\n",
    "\n",
    "print(f\"Total optimal tour length is {OPTIMAL}.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_graph(graph, only_nodes=False):\n",
    "    \"\"\"\n",
    "    Helper method for drawing TSP (tour) graphs.\n",
    "    \"\"\"\n",
    "    fig, ax = plt.subplots(figsize=(12, 6))\n",
    "\n",
    "    if only_nodes:\n",
    "        nx.draw_networkx_nodes(graph, DATA.node_coords, node_size=25, ax=ax)\n",
    "    else:\n",
    "        nx.draw_networkx(\n",
    "            graph, DATA.node_coords, node_size=25, with_labels=False, ax=ax\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_graph(DATA.get_graph(), only_nodes=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use the ALNS meta-heuristic, we need to have destroy and repair operators that work on a proposed solution, and a way to describe such a solution in the first place.\n",
    "Let's start with the solution state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class TspState:\n",
    "    \"\"\"\n",
    "    Solution class for the TSP problem. It has two data members, nodes, and edges.\n",
    "    nodes is a list of IDs. The edges data member, then, is a mapping from each node\n",
    "    to their only outgoing node.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, nodes, edges):\n",
    "        self.nodes = nodes\n",
    "        self.edges = edges\n",
    "\n",
    "    def objective(self):\n",
    "        \"\"\"\n",
    "        The objective function is simply the sum of all individual edge lengths,\n",
    "        using the rounded Euclidean norm.\n",
    "        \"\"\"\n",
    "        return sum(DIST[node, self.edges[node]] for node in self.nodes)\n",
    "\n",
    "    def to_graph(self):\n",
    "        \"\"\"\n",
    "        NetworkX helper method.\n",
    "        \"\"\"\n",
    "        graph = nx.Graph()\n",
    "\n",
    "        for node in self.nodes:\n",
    "            graph.add_node(node, pos=DATA.node_coords[node])\n",
    "\n",
    "        for node_from, node_to in self.edges.items():\n",
    "            graph.add_edge(node_from, node_to)\n",
    "\n",
    "        return graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Destroy operators\n",
    "\n",
    "Destroy operators break parts of a solution down, leaving an incomplete state. This is the first part of each iteration of the ALNS meta-heuristic; the incomplete solution is subsequently repaired by any one repair operator. We will consider three destroy operators: **worst removal**, **path removal** and **random removal**. We will also use a separate parameter, the degree of destruction, to control the extent of the damage done to a solution in each step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DEGREE_OF_DESTRUCTION = 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def edges_to_remove(state):\n",
    "    return int(len(state.edges) * DEGREE_OF_DESTRUCTION)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def worst_removal(current, rng):\n",
    "    \"\"\"\n",
    "    Worst removal iteratively removes the 'worst' edges, that is,\n",
    "    those edges that have the largest distance.\n",
    "    \"\"\"\n",
    "    destroyed = copy.deepcopy(current)\n",
    "\n",
    "    worst_edges = sorted(\n",
    "        destroyed.nodes, key=lambda node: DIST[node, destroyed.edges[node]]\n",
    "    )\n",
    "\n",
    "    for idx in range(edges_to_remove(current)):\n",
    "        del destroyed.edges[worst_edges[-(idx + 1)]]\n",
    "\n",
    "    return destroyed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def path_removal(current, rng):\n",
    "    \"\"\"\n",
    "    Removes an entire consecutive sub-path, that is, a series of\n",
    "    contiguous edges.\n",
    "    \"\"\"\n",
    "    destroyed = copy.deepcopy(current)\n",
    "\n",
    "    node_idx = rng.choice(len(destroyed.nodes))\n",
    "    node = destroyed.nodes[node_idx]\n",
    "\n",
    "    for _ in range(edges_to_remove(current)):\n",
    "        node = destroyed.edges.pop(node)\n",
    "\n",
    "    return destroyed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_removal(current, rng):\n",
    "    \"\"\"\n",
    "    Random removal iteratively removes random edges.\n",
    "    \"\"\"\n",
    "    destroyed = copy.deepcopy(current)\n",
    "\n",
    "    for idx in rng.choice(\n",
    "        len(destroyed.nodes), edges_to_remove(current), replace=False\n",
    "    ):\n",
    "        del destroyed.edges[destroyed.nodes[idx]]\n",
    "\n",
    "    return destroyed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repair operators\n",
    "\n",
    "We implement a simple, **greedy repair** strategy. It determines a set of nodes that are currently not visited, and then links these up to the tour such that it forms one cycle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def would_form_subcycle(from_node, to_node, state):\n",
    "    \"\"\"\n",
    "    Ensures the proposed solution would not result in a cycle smaller\n",
    "    than the entire set of nodes. Notice the offsets: we do not count\n",
    "    the current node under consideration, as it cannot yet be part of\n",
    "    a cycle.\n",
    "    \"\"\"\n",
    "    for step in range(1, len(state.nodes)):\n",
    "        if to_node not in state.edges:\n",
    "            return False\n",
    "\n",
    "        to_node = state.edges[to_node]\n",
    "\n",
    "        if from_node == to_node and step != len(state.nodes) - 1:\n",
    "            return True\n",
    "\n",
    "    return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def greedy_repair(current, rng):\n",
    "    \"\"\"\n",
    "    Greedily repairs a tour, stitching up nodes that are not departed\n",
    "    with those not visited.\n",
    "    \"\"\"\n",
    "    visited = set(current.edges.values())\n",
    "\n",
    "    # This kind of randomness ensures we do not cycle between the same\n",
    "    # destroy and repair steps every time.\n",
    "    shuffled_idcs = rng.permutation(len(current.nodes))\n",
    "    nodes = [current.nodes[idx] for idx in shuffled_idcs]\n",
    "\n",
    "    while len(current.edges) != len(current.nodes):\n",
    "        node = next(node for node in nodes if node not in current.edges)\n",
    "\n",
    "        # Computes all nodes that have not currently been visited,\n",
    "        # that is, those that this node might visit. This should\n",
    "        # not result in a subcycle, as that would violate the TSP\n",
    "        # constraints.\n",
    "        unvisited = {\n",
    "            other\n",
    "            for other in current.nodes\n",
    "            if other != node\n",
    "            if other not in visited\n",
    "            if not would_form_subcycle(node, other, current)\n",
    "        }\n",
    "\n",
    "        # Closest visitable node.\n",
    "        nearest = min(unvisited, key=lambda other: DIST[node, other])\n",
    "\n",
    "        current.edges[node] = nearest\n",
    "        visited.add(nearest)\n",
    "\n",
    "    return current"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial solution\n",
    "\n",
    "We start from a very simple, greedily constructed initial solution. This solution is not good (as can clearly be seen below), but it is feasible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = rnd.default_rng(SEED)\n",
    "state = TspState(CITIES, {})\n",
    "\n",
    "init_sol = greedy_repair(state, rng)\n",
    "print(f\"Initial solution objective is {init_sol.objective()}.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_graph(init_sol.to_graph())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heuristic solution\n",
    "\n",
    "Here we perform the ALNS procedure. The heuristic is given 60 seconds of runtime."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alns = ALNS(rnd.default_rng(SEED))\n",
    "\n",
    "alns.add_destroy_operator(random_removal)\n",
    "alns.add_destroy_operator(path_removal)\n",
    "alns.add_destroy_operator(worst_removal)\n",
    "\n",
    "alns.add_repair_operator(greedy_repair)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "select = RouletteWheel([3, 2, 1, 0.5], 0.8, 3, 1)\n",
    "accept = HillClimbing()\n",
    "stop = MaxRuntime(60)\n",
    "\n",
    "result = alns.iterate(init_sol, select, accept, stop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "solution = result.best_state\n",
    "objective = solution.objective()\n",
    "pct_diff = 100 * (objective - OPTIMAL) / OPTIMAL\n",
    "\n",
    "print(f\"Best heuristic objective is {objective}.\")\n",
    "print(\n",
    "    f\"This is {pct_diff:.1f}% worse than the optimal solution, which is {OPTIMAL}.\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, ax = plt.subplots(figsize=(12, 6))\n",
    "result.plot_objectives(ax=ax, lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at the solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_graph(solution.to_graph())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "In the code above we implemented a very simple heuristic for the TSP, using the ALNS meta-heuristic framework. We did not tinker too much with the various hyperparameters available on the ALNS implementation, but even for these relatively basic heuristic methods and workflow we find a very good result - just a few percent worse than the optimal tour.\n",
    "\n",
    "This notebook showcases how the ALNS library may be put to use to construct powerful, efficient heuristic pipelines from simple, locally greedy operators."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
